How to read this handbook:
This handbook in an HTML file. It is not online, you are only using your web browser to view this local file.
This handbook is best viewed with Google Chrome. Some functions may not work in other browsers.
Use tabs on the right to hide/view code. See ‘Copy to clipboard’ icon in the upper-right of each code section
Version
The latest version of this handbook can be found at this github repository.
Style
dplyr::mutate()). Most times, the code will run without the package being called (without dplyr::), but we do this just to be clear to the reader which package is being used.mutate()).Note types
FOR EXAMPLE: This is a boxed example
TIP: This could also be a note.
NOTE: This is a note
TIP: This is a tip.
CAUTION: This is a cautionary note.
DANGER: This is a warning.
Here the datasets used in this handbook will be described and will be “downloadable” via link (the files will be stored within the HTML, so available offline as well)
Maintainer: Neale Batra (neale.batra@gmail.com)
Code contributors: …
Data contributors: outbreaks package
Content provided by these people…. a…b…c…d…
Review provided by these people…
Some of this material comes from the R4Epis website, which was also made by some of the same people…
RECON
Photo credits (logo)
These tabs will provide basic R skills
How to install R
How to install R Studio
Other things you may need to install:
* TinyTeX
* Pandoc
* RTools
First, open RStudio. As their icons can look very similar, be sure you are opening RStudio and not R.
For RStudio to function you must also have R installed on the computer (see this section for installation instructions).
RStudio is an interface (GUI) for easier use of R. You can think of R as being the engine of a vehicle, doing the crucial work, and RStudio as the body of the vehicle (with seats, accessories, etc.) that helps you actually use the engine to move forward!
By default RStudio displays four rectangle panes.
TIP: If your RStudio displays only one left pane it is because you have no scripts open yet.
The R Console Pane
The R Console, by default the left or lower-left pane in R Studio, is the home of the R “engine”. This is where the commands are actually run and non-graphic outputs and error/warning messages appear. You can directly enter and run commands in the R Console, but realize that these commands are not saved as they are when running commands from a script.
If you are familiar with Stata, the R Console is like the Command Window and also the Results Window.
The Source Pane
This pane, by default in the upper-left, is space to edit and run your scripts. This pane can also display datasets (data frames) for viewing.
For Stata users, this pane is similar to your Do-file and Data Editor windows.
The Environment Pane
This pane, by default the upper-right, is most often used to see brief summaries of objects in the R Environment in the current session. These objects could include imported, modified, or created datasets, parameters you have defined (e.g. a specific epi week for the analysis), or vectors or lists you have defined during analysis (e.g. names of regions). Click on the arrow next to a dataframe name to see its variables.
In Stata, this is most similar to Variables Manager window.
Plots, Packages, and Help Pane
The lower-right pane includes several tabs including plots (display of graphics including maps), help, a file library, and available R packages (including installation/update options).
This pane contains the Stata equivalents of the Plots Manager and Project Manager windows.
Change RStudio settings and appearance in the Tools drop-down menu, by selecting Global Options
R scripts (vs. typing in the console)
* Advantages (reproducability) * General sequence (into, load packages, load data, clean data, conduct analysis, save results) * Commenting
These tabs cover how to use R working directories, and how this changes when you are working within an R project. The working directory is the root file location used by R for your work.
By default, it will save new files and outputs to this location, and will look for files to import (e.g. datasets) here as well.
NOTE: If using an [R project](#rproject), the working directory will default to the R project root folder **IF** you open RStudio by clicking open the R project (the file with .rproj extension))
Use the command setwd() with the filepath in quotations, for example: setwd("C:/Documents/R Files")
CAUTION: If using an RMarkdown script be aware of the following:
In an R Markdown script, the default working directory is the folder the Rmarkdown file (.Rmd) is saved to. If you want to change this, you can use setwd() as above, but know the change will only apply to that specific code chunk.
To change the working directory for all code chunks in an R markdown, edit the setup chunk to add the root.dir = parameter, such as below:
Setting your working directory manually (point-and-click)
From RStudio click: Session / Set Working Directory / Choose Directory (you will have to do this each time you open RStudio)
How things change in an R project
Everything in R is an object. These sections will explain:
<-)Everything you store in R - datasets, variables, a list of village names, a total population number, even outputs such as graphs - are objects which are assigned a name and can be referenced in later commands.
An object exists when you have assigned it a value (see the assignment section below). When it is assigned a value, the object appears in the Environment (see the upper right pane of RStudio). It can then be operated upon, manipulated, changed, and re-defined.
<-)Create objects by assigning them a value with the <- operator.
You can think of the assignment operator <- as the words “is defined as”. Assignment commands generally follow a standard order:
object_name <- value (or process/calculation that produce a value)
EXAMPLE: You may want to record the current epidemiological reporting week as an object for reference in later code. In this example, the object
reporting_weekis created when it is assigned the character value"2018-W10"(the quote marks make these a character value).
The objectreporting_weekwill then appear in the RStudio Environment pane (upper-right) and can be referenced in later commands.
See the R commands and their output in the boxes below.
reporting_week <- "2018-W10" # this command creates the object reporting_week by assigning it a value
reporting_week # this command prints the current value of reporting_week object in the console
## [1] "2018-W10"NOTE: Note the [1] in the R console output is simply indicating that you are viewing the first item of the output
CAUTION: An object’s value can be over-written at any time by running an assignment command to re-define its value. Thus, the order of the commands run is very important.
The following command will re-define the value of reporting_week:
reporting_week <- "2018-W51" # assigns a NEW value to the object reporting_week
reporting_week # prints the current value of reporting_week in the console
## [1] "2018-W51"Datasets are also objects and must be assigned names when they are imported.
In the code below, the object linelist_raw is created and assigned the value of a CSV file imported with the rio package.
# linelist_raw is created and assigned the value of the imported CSV file
linelist <- rio::import("my_linelist.csv")You can read more about importing and exporting datasets with the section on importing data.
CAUTION: A quick note on naming of objects:
Objects can be a single piece of data (e.g. my_number <- 24), or they can consist of structured data.
The graphic below, sourced from this online R tutorial shows some common data structures and their names. Not included in this image is spatial data, which is discussed in the GIS section.
In epidemiology (and particularly field epidemiology), you will most commonly encounter data frames and vectors:
| Common structure | Explanation | Example from templates |
|---|---|---|
| Vectors | A container for a sequence of singular objects, all of the same class (e.g. numeric, character). | “Variables” (columns) in data frames are vectors (e.g. the variable age_years). |
| Data Frames | Vectors (e.g. columns) that are bound together that all have the same number of rows. | linelist_raw and linelist_cleaned are both data frames. |
Note that to create a vector that “stands alone”, or is not part of a data frame (such as a list of location names), the function c() is often used:
list_of_names <- c("Ruhengeri", "Gisenyi", "Kigali", "Butare")
$)Vectors within a data frame (variables in a dataset) can be called, referenced, or created using the $ symbol. The $ symbol connects the name of the column to the name of its data frame. The $ symbol must be used, otherwise R will not know where to look for or create the column.
# Retrieve the length of the vector age_years
length(linelist$age) # (age is a variable in the linelist data frame)By typing the name of the data frame followed by $ you will also see a list of all variables in the data frame. You can scroll through them using your arrow key, select one with your Enter key, and avoid spelling mistakes!
All the objects stored in R have a class which tells R how to handle the object. There are many possible classes, but common ones include:
| Class | Explanation | Examples |
|---|---|---|
| Character | These are text/words/sentences “within quotation marks”. Math cannot be done on these objects. | “Character objects are in quotation marks” |
| Numeric | These are numbers and can include decimals. If within quotation marks the will be considered character. | 23.1 or 14 |
| Integer | Numbers that are whole only (no decimals) | -5, 14, or 2000 |
| Factor | These are vectors that have a specified order or hierarchy of values | Variable msf_involvement with ordered values N, S, SUB, and U. |
| Date | Once R is told that certain data are Dates, these data can be manipulated and displayed in special ways. See the page on Dates for more information. | 2018-04-12 or 15/3/1954 or Wed 4 Jan 1980 |
| Logical | Values must be one of the two special values TRUE or FALSE (note these are not “TRUE” and “FALSE” in quotation marks) | TRUE or FALSE |
| data.frame | A data frame is how R stores a typical dataset. It consists of vectors (columns) of data bound together, that all have the same number of observations (rows). | The example AJS dataset named linelist_raw contains 68 variables with 300 observations (rows) each. |
You can test the class of an object by feeding it to the function class(). Note: you can reference a specific column within a dataset using the $ notation to separate the name of the dataset and the name of the column.
class(linelist$age) # class should be numeric
## [1] "numeric"
class(linelist$gender) # class should be character
## [1] "character"Often, you will need to convert objects or variables to another class.
| Function | Action |
|---|---|
as.character() |
Converts to character class |
as.numeric() |
Converts to numeric class |
as.integer() |
Converts to integer class |
as.Date() |
Converts to Date class - Note: see section on dates for details |
as.factor() |
Converts to factor - Note: re-defining order of value levels requires extra arguments |
Here is more online material on classes and data structures in R.
This section on functions explains:
* What a function is and how they work
* What arguments are
* What packages are
* How to get help understanding a function
A function is like a machine that receives inputs, does some action with those inputs, and produces an output.
What the output is depends on the function.
Functions typically operate upon some object placed within the function’s parentheses. For example, the function sqrt() calculates the square root of a number:
Functions can also be applied to variables in a dataset. For example, when the function summary() is applied to the numeric variable age in the dataset linelist (what’s the $ symbol?), the output is a summary of the variable’s numeric and missing values.
summary(linelist$age)
## Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
## 2.00 47.25 60.00 56.91 72.00 91.00 2NOTE: Behind the scenes, a function represents complex additional code that has been wrapped up for the user into one easy command.
Functions often ask for several inputs, called arguments, located within the parentheses of the function, usually separated by commas.
For example, this
age_pyramid()command produces an age pyramid graphic based on defined age groups and a binary split variable, such asgender. The function is given three arguments within the parentheses, separated by commas. The values supplied to the arguments establishlinelistas the data frame to use,age_groupas the variable to count, andgenderas the binary variable to use for splitting the pyramid by color.
NOTE: For this example, in the background we have created a new variable called “age_group”. To learn how to create new variable see that section of this handbook
# Creates an age pyramid by specifying the dataframe, age group variable, and a variable to split the pyramid
apyramid::age_pyramid(data = linelist, age_group = "age_group", split_by = "gender")The first half of an argument assignment (e.g.
data =) does not need to be specified if the arguments are written in a specific order (specified in the function’s documentation). The below code produces the exact same pyramid as above, because the function expects the argument order: data frame,age_groupvariable,split_byvariable.
# This command will produce the exact same graphic as above
apyramid::age_pyramid(linelist, "age_group", "gender")A more complex age_pyramid() command might include the optional arguments to:
proportional = TRUE when the default is FALSE)pal = is short for “palette” and is supplied with a vector of two color names. See the objects page for how the function c() makes a vector)NOTE: For arguments specified with an equals symbol (e.g. coltotals = ...), their order among the arguments is not important (must still be within the parentheses and separated by commas).
Packages contain functions.
On installation, R contains “base” functions that perform common elementary tasks. But many R users create specialized functions, which are verified by the R community and which you can download as a package for your own use.
One of the more challenging aspects of R is that there are often many functions or packages to choose from to complete a given task.
Functions are contained within packages which can be downloaded (“installed”) to your computer from the internet. Once a package is downloaded, you access its functions by loading the package with the library() command at the beginning of each R session.
NOTE: While you only have to install a package once, you must load it at the beginning of every R session using library() command, or an alternative like pacman’s p_load() function.
Think of R as your personal library: When you download a package your library gains a book of functions, but each time you want to use a function in that book, you must borrow that book from your library.
For clarity in this handbook, functions are usually preceeded by the name of their package using the :: symbol in the following way:
package_name::function_name()
Once a package is loaded for a session, this explicit style is not necessary. One can just use function_name(). However giving the package name is useful when a function name is common and may exist in multiple packages (e.g. plot()).
Using the package name will also load the package if it is not already loaded.
# This command uses the package "rio" and its function "import()" to import a dataset
linelist <- rio::import("linelist.xlsx", which = "Sheet1")Dependencies
Packages often depend on other packages, and these are called “dependencies”. When a package is installed from CRAN, it will typically also install its dependenices.
To read more about a function, you can try searching online for resources OR search in the Help tab of the lower-right RStudio pane.
%>%)Two general approaches to R coding are:
Simply explained, the pipe operator (%>%) passes an intermediate output from one function to the next.
You can think of it as saying “then”. Many functions can be linked together with %>%.
Piping emphasizes a sequence of actions, not the object the actions are being performed on
Best when a sequence of actions must be performed on one object
from magrittr. Included in dplyr and tidyverse
Makes code more clean and easier to read, intuitive
express a sequence of operations
the object is altered and then passed on to the next function
Example:
# A fake example of how to bake a care using piping syntax
cake <- flour %>% # to define cake, start with flour, and then...
left_join(eggs) %>% # add eggs
left_join(oil) %>% # add oil
left_join(water) %>% # add water
mix_together(utensil = spoon, minutes = 2) %>% # mix together
bake(degrees = 350, system = "fahrenheit", minutes = 35) %>% # bake
let_cool() # let it cool downhttps://cfss.uchicago.edu/notes/pipes/#:~:text=Pipes%20are%20an%20extremely%20useful,code%20and%20combine%20multiple%20operations.
Piping is not a base function. To use piping, the dplyr package must be installed and loaded. Near the top of every template script is a code chunk that installs and loads the necessary packages, including dplyr. You can read more about piping in the documentation.
CAUTION: Remember that even when using piping to link functions, if the assignment operator (<-) is present, the object to the left will still be over-written (re-defined) by the right side.
Better if:
* You need to manipulate multiple objects
* There are intermediate steps that are meaningful and deserve separate object names
as changes are made - still handy to know
Risks: creating new objects for each step - lots of objects. If you use the wrong one you might not know. naming can be confusing, errors not easily detectable
either name each intermediate object, or overwrite the original, or combine all the functions together. all come with risks
https://style.tidyverse.org/pipes.html
# a fake example of how to bake a cake using this method (defining intermediate objects)
batter_1 <- left_join(flour, eggs)
batter_2 <- left_join(batter_1, oil)
batter_3 <- left_join(batter_2, water)
batter_4 <- mix_together(object = batter_3, utensil = spoon, minutes = 2)
cake <- bake(batter_4, degrees = 350, system = "fahrenheit", minutes = 35)
cake <- let_cool(cake)Combine all functions together - also difficult to read
# an example of combining/nesting mutliple functions together - difficult to read
cake <- let_cool(bake(mix_together(batter_3, utensil = spoon, minutes = 2), degrees = 350, system = "fahrenheit", minutes = 35))This section details operators in R, such as:
* Relational operators (less than, equal too..)
* Logical operators (and, or…)
* Missingness
* Mathematical operators (+, -, /…) * The %in% operator
Relational operators compare values and are often used when defining new variables and subsets of datasets. Here are the common relational operators in R:
| Function | Operator | Example | Example Result |
|---|---|---|---|
| Equal to | == |
"A" == "a" |
FALSE (because R is case sensitive) Note that == (double equals) is different from = (single equals), which acts like the assignment operator <- |
| Not equal to | != |
2 != 0 |
TRUE |
| Greater than | > |
4 > 2 |
TRUE |
| Less than | < |
4 < 2 |
FALSE |
| Greater than or equal to | >= |
6 >= 4 |
TRUE |
| Less than or equal to | <= |
6 <= 4 |
FALSE |
| Value is missing | is.na() |
is.na(7) |
FALSE (see section on missing values) |
| Value is not missing | !is.na() |
!is.na(7) |
TRUE |
Logical operators, such as AND and OR, are often used to connect relational operators and create more complicated criteria. Complex statements might require parentheses ( ) for grouping and order of application.
| Function | Operator |
|---|---|
| AND | & |
| OR | | (vertical bar) |
| Parentheses | ( ) Used to group criteria together and clarify order |
For example, below, we have a linelist with two variables we want to use to create our case definition, hep_e_rdt, a test result and other_cases_in_hh, which will tell us if there are other cases in the household. The command below uses the function case_when() to create the new variable case_def such that:
linelist_cleaned <- linelist_cleaned %>%
mutate(case_def = case_when(
is.na(hep_e_rdt) & is.na(other_cases_in_hh) ~ NA_character_,
hep_e_rdt == "Positive" ~ "Confirmed",
hep_e_rdt != "Positive" & other_cases_in_hh == "Yes" ~ "Probable",
TRUE ~ "Suspected"
))| Criteria in example above | Resulting value in new variable “case_def” |
|---|---|
If the value for variables hep_e_rdt and other_cases_in_hh are missing |
NA (missing) |
If the value in hep_e_rdt is “Positive” |
“Confirmed” |
If the value in hep_e_rdt is NOT “Positive” AND the value in other_cases_in_hh is “Yes” |
“Probable” |
| If one of the above criteria are not met | “Suspected” |
{{% notice tip %}} Note that R is case-sensitive, so “Positive” is different than “positive”… {{% /notice %}}
In R, missing values are represented by the special value NA (capital letters N and A - not in quotation marks). If you import data that records missing data in another way (e.g. 99, “Missing”, or .), you may want to re-code those values to NA.
To test whether a value is NA, use the special function is.na(), which returns TRUE or FALSE.
Mathematical operators are often used to perform addition, division, to create new columns, etc. Below are common mathematical operators in R. Whether you put spaces around the operators is not important.
| Objective | Example in R |
|---|---|
| addition | 2 + 3 |
| subtraction | 2 - 3 |
| multiplication | 2 * 3 |
| division | 30 / 5 |
| exponent | 2^3 |
| order of operations | ( ) |
%in%This section describes the several ways to install a package:
* Via the online package repository (CRAN)
* From a ZIP file
* From Github
This section explains:
* General syntax for writing R code
* Code assists
* the difference between errors and warnings
Common errors and warnings and their solutions can be found in X section (TODO).
A few things to remember when writing commands in R, to avoid errors and warnings:
Variable_A is different from variable_AAny script (RMarkdown or otherwise) will give clues when you have made a mistake. For example, if you forgot to write a comma where it is needed, or to close a parentheses, RStudio will raise a flag on that line, on the right side of the script, to warn you.
(/images/Warnings_and_Errors.png)
When a command is run, the R Console may show you warning or error messages in red text.
A warning means that R has completed your command, but had to take additional steps or produced unusual output that you should be aware of.
An error means that R was not able to complete your command.
Look for clues:
The error/warning message will often include a line number for the problem.
If an object “is unknown” or “not found”, perhaps you spelled it incorrectly, forgot to call a package with library(), or forgot to re-run your script after making changes.
If all else fails, copy the error message into Google along with some key terms - chances are that someone else has worked through this already!
Introduction to importing data
The key package we recommend for importing data is: rio. The package rio offers the useful function import() that can import many types of files into R.
The alternative to using rio would be to use functions from several other packages that are specific to a type of file (e.g. read.csv(), read.xlsx(), etc.). While these alternatives can be difficult to remember, always using rio::import() is relatively easy.
Optionally, the package here can be used in conjunction with rio. It locates files on your computer via relative pathways, usually within the context of an R project. Relative pathways are relative from a designated folder location, so that pathways listed in R code will not break when the script is run on a different computer.
This code chunk shows the loading of packages for importing data.
import()When you import a dataset, you are doing the following:
[TODO importing image]
The function import() from the package rio makes it easy to import many types of data files.
import() uses the file’s extension (e.g. .xlsx, .csv, .dta, etc.) to appropriately import the file. Any optional arguments specific to the filetype can be supplied as well.
You can read more about the rio package in this online vignette
Manually typing a filepath is prone to error, but is sometime the best way (for example, if referencing files from a shared/network drive). The function import() (from the package rio) accepts a filepath in quotes.
A few things to note:
If importing a specific sheet from an Excel file, include the sheet name in the which = argument of import(). For example:
# A demonstration showing how to import a specific Excel sheet
my_data <- rio::import("my_excel_file.xlsx", which = "Sheetname")If using the here() method to provide a relative pathway to import(), you can still indicate a specific sheet by adding the which = argument after the closing parenthese of the here() function.
You can import data manually by either:
file.choose() (leaving the parentheses empty) to trigger appearance of a pop-up window that allows the user to manually select the file from their computer. For example:here())Relative filepaths differ from static filepaths in that they are relative from a R project root directory. For example:
import("C:/Users/nsbatra/My Documents/R files/epiproject/data/linelists/ebola_linelist.xlsx")
import(here("data", "linelists", "ebola_linelist.xlsx"))
The package here and it’s function here() facilitate relative pathways.
here() works best within R projects. When the here package is first loaded (library(here)), it automatically considers the top-level folder of your R project as “here” - a benchmark for all other files in the project.
Thus, in your script, if you want to import or reference a file saved in your R project’s folders, you use the function here() to tell R where the file is in relation to that benchmark.
If you are unsure where “here” is set to, run the function here() with the empty brackets:
Below is an example of importing the file “fluH7N9_China_2013.csv” which is located in the benchmark “here” folder. All you have to do is provide the name of the file in quotes (with the appropriate ending).
If the file is within a subfolder - let’s say a “data” folder - write these folder names in quotes, separated by commas, as below:
Using the here() command produces a character filepath, which can then processed by the import() function.
# the filepath
here("data", "fluH7N9_China_2013.csv")
# the filepath is given to the import() function
linelist <- import(here("data", "fluH7N9_China_2013.csv"))NOTE: You can still import a specific sheet of an excel file as noted elsewhere. The here() command is only supplying the filepath.
By default, R expect manual entry of data in vectors (columns) and combining them into a data frame: NOTE: all vectors must be the same length!
# define each vector (column) separately, each with their own name
PatientID <- c(235, 452, 778, 111)
Treatment <- c("Yes", "No", "Yes", "Yes")
Death <- c(1, 0, 1, 0)
# combine the columns into one data frame, by referencing the vector names
manual_entry_cols <- data.frame(PatientID, Treatment, Death)And now we display the new dataset:
Use the tribble function from the tibble package from the tidverse.
https://tibble.tidyverse.org/reference/tribble.html
Note that each column must contain only one kind of data (character, numeric, etc.). You can create a list column, however (see documentation) You can use tabs and spacing to make the data entry more intuitive and readable. For example:
# create the dataset manually by row
manual_entry_rows <- tibble::tribble(
~colA, ~colB,
"a", 1,
"b", 2,
"c", 3
)And now we display the new dataset:
x <- read.table(file = “clipboard”, sep = “t”, header=TRUE) x
The following packages are recommended for working with dates:
as.Date()The standard, base R function to convert an object or variable to class Date is as.Date() (not capitalization).
as.Date() requires that the user specify the existing format so that it is able to convert and store each element correctly (day, month, year, etc.).
If used on a variable, as.Date() therefore requires that all the character date values be in the same format before converting. Read more about using as.Date() here.
It can be easiest to first convert the variable to character class, and then convert to date class:
as.character()as.Date()Within the as.Date() function, you must use the format= argument to tell R which characters are which date components - which characters refer to the month, the day, and the year. If your values are already in one of R’s standard formats (YYYY-MM-DD or YYYY/MM/DD) the format= argument is not necessary.
For example, if your character dates are in the format DD/MM/YYYY, like “24/04/1968”, then your command to turn the values into dates will be as below. Putting the format in quotation marks is necessary.
TIP: The format= argument is not telling R the format you want the dates to be, but rather how to identify the date parts as they are before you run the command.
TIP: Also, be sure that in the format= argument you use the date-part separator (e.g. /, -, or space) that is present in your dates.
The as.character() and as.Date() commands can be wrapped together such as:
linelist_cleaned$date_of_onset <- as.Date(as.character(linelist_cleaned$date_of_onset), format = "%d/%m/%Y")or within piping (see REF):
# Convert the variable to class Date by providing the format of the variable
linelist_cleaned$date_of_onset <- as.Date(linelist_cleaned$date_of_onset, format="%Y-%m-%d")
# Check the class of the variable again
class(linelist_cleaned$date_of_onset) Once the values are in class Date, R will present them in it’s standard format, which is YYYY-MM-DD.
guess_dates()The function guess_dates() attempts to read a “messy” date variable containing dates in many different formats and convert the dates to a standard format. You can read more about guess_dates(), which is in the linelist package.
For example:
guess_dateswould see the following dates “03 Jan 2018”, “07/03/1982”, and “08/20/85” and convert them in the class Date to: 2018-01-03, 1982-03-07, and 1985-08-20.
#guess_dates(c("03 Jan 2018", "07/03/1982", "08/20/85")) # guess_dates() not yet available for R 4.0.2Some optional arguments for guess_dates() that you might include are:
error_tolerance - The proportion of entries which cannot be identified as dates to be tolerated (defaults to 0.1 or 10%)last_date - the last valid date (defaults to current date)first_date - the first valid date. Defaults to fifty years before the last_date.# An example using guess_dates on the variable dtdeath
linelist_cleaned$dtdeath <- linelist::guess_dates(linelist_cleaned$dtdeath)
# An example from the template using guess_dates over multiple date variables, with piping, error tolerance of 50%, and the earliest accepted date of 1 Jan 2016.
linelist_cleaned <- linelist_cleaned %>%
mutate_at(vars(matches("date|Date")), linelist::guess_dates,
error_tolerance = 0.5, first_date = "2016-01-01")Excel stores dates as the number of days since December 30, 1899. If the dataset you imported from Excel has a date variable showing numbers like “41369”… use the as.Date() function to convert, but instead of supplying a format as above, supply an origin date. Note that the origin date must be given in the default date format for R (“YYYY-MM-DD”).
Once dates are the correct class, you often want them to display differently. For example, “Monday 05 Jan” instead of 2018-01-05. You can do this with the function format() in a similar was to as.Date(). Read more about it in this tutorial
The difference between dates can be calculated by:
The templates use the very flexible package aweek to set epidemiological weeks. You can read more about it on the RECON website
See the section on epicurves.
Text about cleaning data, approaches, etc. renaming replace missing with dealing with cases (all lower, etc) case_when() factors
Variable names are used so often, it is best to have them be “clean” (no spaces, no unusual characters, etc.)
The function clean_names() from the package janitor is very useful.
https://cran.r-project.org/web/packages/janitor/vignettes/janitor.html#cleaning
base R method
Using dplyr mutate()
example
case_when())For example, creating age groups cut()
case_when()
age_categories() (R4Epis package)
Time zones Shifting time (adding hours) example of WHO receiving reports of COVID from different time zones
Missing if… na_if()
Replace
Santa Clara County example - COVID contact tracing data - classification of multiple phone call records from same person into the highest category. (classify all as the highest of the group)
dealing with missing data percent missing over time etc.
Or change in percent of anything (X) over time, really.
lines <- linelist %>%
mutate(date_of_onset = as.Date(date_of_onset, format = "%d/%m/%Y"),
week = aweek::week2date(aweek::date2week(date_of_onset))) %>%
group_by(week) %>%
summarize(n_obs = n(),
dt_hosp_missing = sum(date_of_hospitalisation == "" | is.na(date_of_hospitalisation)),
dt_hosp_p_miss = dt_hosp_missing / n_obs,
outcome_missing = sum(outcome == "" | is.na(outcome)),
outcome_p_miss = outcome_missing / n_obs) %>%
reshape2::melt(id.vars = c("week")) %>%
filter(grepl("_p_", variable)) %>%
ggplot()+
geom_line(aes(x = week, y = value, group = variable, color = variable), size = 1, stat = "identity")+
labs(title = "Missingness in variables, as proportion of ",
#subtitle = str_glue("As of {format(report_date, '%d %b')}"),
x = "Week",
y = "Proportion missing",
fill = "CalREDIE Variable") +
scale_color_discrete(name = "Variable", labels = c("Date of Hospitalization Missing", "Outcome Missing"))+
scale_y_continuous(breaks = c(seq(0,1,0.1)))
#theme_cowplot()#+
#theme(legend.position = element_text("none"))
lines(pivoting/melting etc.) Transforming datasets from wide-to-long, or long-to-wide…
Transforming a dataset from wide to long
We start with data that is in a wide format, e.g. our linelist.
pivot_longer()dplyr pivot_wider()
Tidyverse - grouping by values
.drop=F in group_by() command
group_by()aggregate()This section includes:
* Basic descriptive statistics on numeric variables
* …
* …
Note the argument na.rm=T, which removes missing values from the calculation.
If missing values are not excluded, the returned value will be NA (missing).
Note the argument na.rm=T, which removes missing values from the calculation.
If missing values are not excluded, the returned value will be NA (missing).
Note the argument na.rm=T, which removes missing values from the calculation.
If missing values are not excluded, the returned value will be NA (missing).
Note the argument na.rm=T, which removes missing values from the calculation.
If missing values are not excluded, the returned value will be NA (missing).
Note the argument na.rm=T, which removes missing values from the calculation.
If missing values are not excluded, the returned value will be NA (missing).
Frequency table of 1 and 2 categorical variables
table(linelist$province)
##
## Anhui Beijing Fujian Guangdong Hebei Henan Hunan Jiangsu
## 4 3 5 1 1 4 2 28
## Jiangxi Shandong Shanghai Taiwan Zhejiang
## 6 2 33 1 46
x <- table(linelist$province, linelist$gender)
#janitor::adorn_totals(x)A table with 3 variables
table_3vars <- table(linelist$province, linelist$gender, linelist$outcome)
ftable(table_3vars)
## Death Recover
##
## Anhui f 1 0
## m 1 1
## Beijing f 0 1
## m 0 1
## Fujian f 0 0
## m 0 3
## Guangdong f 0 0
## m 0 0
## Hebei f 1 0
## m 0 0
## Henan f 0 0
## m 1 3
## Hunan f 1 0
## m 0 1
## Jiangsu f 2 3
## m 2 7
## Jiangxi f 0 2
## m 1 1
## Shandong f 0 0
## m 0 2
## Shanghai f 3 3
## m 12 10
## Taiwan f 0 0
## m 0 0
## Zhejiang f 1 3
## m 5 5count_data %>%
group_by(District) %>%
summarise(n_obs = n(), # number of observations
range_date = max(Date, na.rm=T) - min(Date, na.rm=T)
)
## # A tibble: 4 x 3
## District n_obs range_date
## <chr> <int> <drtn>
## 1 Basraj 2676 188 days
## 2 Dingo 2073 188 days
## 3 Gargona 2362 188 days
## 4 Nibari 2517 188 daysThese can be created with:
boxplot() function from the graphics package (installed automatically with base R), orggplot() function from the ggplot2 package# with the boxplot() function
graphics::boxplot(age ~ outcome*gender,
data = linelist,
col = c("gold", "darkgreen"),
main = "Made using boxplot()")
# with ggplot2
ggplot(data = linelist %>% filter(!is.na(gender)),
aes(y = age, x = outcome, fill = outcome))+
geom_boxplot()+
ggtitle("Made using ggplot()")+
facet_wrap(~gender)boxplot()Some options with boxplot() shown below are:
* boxplots by group (color specification optional) * violin plots
* boxplot width proportional to sample size
* Horizontal
# boxplot of one numeric variable
boxplot(linelist$age, # numeric variable
main="boxplot", # main title
xlab="Suppliment and Dose") # x-axis label
# by group (formula style)
boxplot(age ~ gender, data=linelist, notch=TRUE, main="boxplot", xlab="Suppliment and Dose")You can have multiple levels of group (e.g. age by outcome AND gender)
Notched “violin plots” are possible. The notch represents the median and X around it (TODO)
# By subgroup (age by outcome AND gender)
boxplot(age ~ outcome * gender,
data=linelist,
col=c("gold","darkgreen"), # colors, in a vector
main="Boxplot by Outcome and Gender", # main title
xlab="Suppliment and Dose") # x-axis label
# Notched (violin plot), and varying width
boxplot(age ~ outcome * gender,
data=linelist,
notch=TRUE, # notch at median
varwidth = TRUE, # width varying by sample size
col=(c("gold","darkgreen")),
main="Notched boxplot, width varying by sample size",
xlab="Suppliment and Dose")
# Horizontal
boxplot(age~outcome,
data=linelist,
horizontal=TRUE, # flip to horizontal
col=(c("gold","darkgreen")),
main="Horizontal boxplot",
xlab="Suppliment and Dose")ggplot()Some options with ggplot() shown below are:
* boxplots by group (color specification optional) * violin plots
* boxplot width proportional to sample size
* Horizontal
# Simple boxplot of one numeric variable
ggplot(data = linelist, aes(y = age))+ # only y variable given (no x variable)
geom_boxplot()+
ggtitle("Simple ggplot() boxplot")
# By group
ggplot(data = linelist, aes(y = age, # numeric variable
x = outcome, # group variable
fill = outcome))+ # fill variable (color of boxes)
geom_boxplot()+ # create the boxplot
ggtitle("ggplot() boxplot by gender") # main title
# Removing missing values, and add color
ggplot(data = linelist %>% filter(!is.na(outcome)), # dataset piped through a filter to retain rows where gender is not missing
aes(y = age, x = outcome, fill= outcome))+ # boxes filled according to gender value
geom_boxplot()+
ggtitle("ggplot() boxplot by gender (missing excluded)")To examine by subgroups, use facet_wrap() (for more see section on ggplot tips).
# By subgroup
ggplot(data = linelist %>% filter(!is.na(gender)), # dataset piped through a filter to retain rows where gender is not missing
aes(y = age, x = outcome, fill=outcome))+
geom_boxplot()+
ggtitle("A ggplot() boxplot")+
facet_wrap(~gender)“Violin plots” can be made simply or very complex:
# Vertical violin plot
ggplot(linelist, aes(x=age, y=outcome, fill = outcome)) +
geom_violin(trim=FALSE)+
coord_flip()
# Add jittered points
ggplot(linelist, aes(x=age, y=outcome, fill = outcome)) +
geom_violin(trim=FALSE)+
coord_flip()+
geom_jitter(shape=16, # points
position=position_jitter(0.2)) # jitter permissible to avoid point overlapquant/quant, quant/cat, cat/cat t-tests odds ratios, mantel-haensel, etc.
Below are the Preparation tabs:
The data preparation involves the following steps, detailed in the following tabs:
Load packages: installs and load the packages required for the scripts
Load data: imports datasets
Clean data: this section contains ad hoc data cleaning, i.e. which is not used in other reports (otherwise cleaning should be done in a dedicated report); this section is also used to create new variables used in the analyses
This code chunk shows the loading of packages required for the analyses.
# Create vector of names of required packages:
packages_epicurve <- c("rio", # File import
"here", # File locator
"tidyverse", # data manipulation
"ggplot2", # Produce plots and graphs
"aweek", # working with dates
"lubridate", # Manipulate dates
"incidence", # an option for epicurves of linelist data
"stringr", # Search and manipulate character strings
"forcats", # working with factors
"RColorBrewer", # Color palettes from colorbrewer2.org
"DT" # produce tables for this html handbook
) ### close vector of packages
# Checks if package is installed, installs if necessary, and loads package for current session
pacman::p_load(packages_epicurve, character.only=TRUE)Two example datasets are used in this section:
If viewing in Google Chrome, you can access these datasets in Microsoft Excel by clicking HERE and HERE. TODO.
The datasets are imported using the import() function from the rio package. See the page on importing data for various ways to import data. Each data set is displayed below as a table for viewing.
For most of this document, the linelist dataset will be used. The aggregated counts dataset will be used at the end.
Review the two datasets and notice the differences
Linelist dataset
# display the linelist data as a table
DT::datatable(linelist, rownames = FALSE, filter="top", options = list(pageLength = 5, scrollX=T) )Aggregated counts dataset
You may want to set certain parameters for production of a report, such as the date for which the data is current (the “data date”). In this case, we set this date as 27 July 2013.
Now we can reference the object data_date into the code and have it reference that date.
Dates are essential to epidemiological curves. You must ensure that R knows which variables contain dates and that they are recognized correctly.
Optionally, it can be nice to identify all the date variables and store their names in a vector. This can be done by individually naming them, or by searching for them by looking for keywords.
Method 1
# Method 1. Define date variables explicitly in a vector
DateVars <- c("date_of_onset",
"date_of_hospitalisation",
"date_of_outcome"
)
DateVars
## [1] "date_of_onset" "date_of_hospitalisation"
## [3] "date_of_outcome"Method 2
# Method 2: Search for date columns
DateVars <- as.character(tidyselect::vars_select(names(linelist), matches("date|Date|dt")))
DateVars
## [1] "date_of_onset" "date_of_hospitalisation"
## [3] "date_of_outcome"
#Note: other search tool options within vars_select include contains() ends_with(), starts_with(), or one_of()Convert the date variables to class “date”. There are a few options that use different packages. Each is explained in the code chunk.
# Method 1. Manually convert each variable (allows flexibility in format)
linelist$date_of_onset <- as.Date(linelist$date_of_onset, format = "%m/%d/%Y")
linelist$date_of_hospitalisation <- as.Date(linelist$date_of_hospitalisation, format = "%m/%d/%Y")
linelist$date_of_outcome <- as.Date(linelist$date_of_outcome, format = "%m/%d/%Y")
# Method 2. Use clean_dates() on the dataset to clean date variables
#linelist <- linelist %>%
# linelist::clean_dates(error_tolerance = 0.1)
# Note: Use guess_dates() function from the linelist package addresses messy dates (different formats within a variable)Verify that each variable was successfully converted to date class by printing statistics and a quick histogram for each one.
# To verify successful conversion of date variables
# Creates list of column numbers of date variables
varNums <- c()
for (varName in DateVars) {
varNum <- match(varName, names(linelist))
varNums <- c(varNums, varNum)
}
# Produce output for each date variable converted
for (varNum in varNums) {
varName <- names(linelist)[varNum] # get name of variable
class <- class(linelist[, varNum]) # get class
missing <- sum(is.na(linelist[, varNum])) # get number missing values
hist(linelist[, varNum], # histogram
breaks = 50,
main = paste0("Histogram of: ", varName, ", Class: ", class, ", Missing: ", missing),
xlab = varName)
}Below are tabs on using the “incidence” package
This section shows variations on the epicurve using the incidence package
These are simple epicurves using the incidence package. The epicurve is assigned to the object “epicurve”, which is then plotted. Remember that incidence::plot() is different to base::plot()
The interval defines how the observations are grouped. Options are all those in the package aweek, including but are not limited to:
* “Monday week” * “2 Monday weeks” * “Sunday week”
* “MMWRweek” (starts on Sunday)
* “Month”
* “Quarter”
* “Year”
First date and last date can also be specified.
# incidence object is created, with data aggregated at one day intervals
epicurve_daily <- incidence::incidence(linelist$date_of_onset, interval = "day")
# If weekly, you can specific the start day
epicurve_weekly <- incidence::incidence(linelist$date_of_onset, interval = "Monday week")
epicurve_3weekly <- incidence::incidence(linelist$date_of_onset, interval = "3 weeks")
# Monthly
epicurve_monthly <- incidence::incidence(linelist$date_of_onset, interval = "month")
# Plot the incidence object
plot(epicurve_daily)
plot(epicurve_weekly)
plot(epicurve_3weekly)
plot(epicurve_monthly)Behind the scenes, incidence is using ggplot(), so you can add aesthetic themes and other lines using the ggplot syntax.
# Set theme elements using ggplot syntax
epicurve_theme <- ggplot2::theme(
axis.text.x = element_text(angle = 45, hjust = 1, vjust = 1),
legend.title = element_blank(),
panel.grid.major.x = element_line(color = "grey60", linetype = 3),
panel.grid.major.y = element_line(color = "grey60", linetype = 3)
)
# Sets labels using ggplot syntax
epicurve_labels <- labs(x = "Week",
y = "Cases (n)",
title = "H5N7 cases by week of onset",
caption = paste0("Source: Linelist data from: ", data_date, "; ", missing_onset, " are missing date of onset and not shown."))
# plot the epicurve with aesthetics
nice_plot <- plot(epicurve_weekly, show_cases = TRUE, border = "black", n_breaks = nrow(epicurve_weekly)) +
scale_y_continuous(expand = c(0, 0)) + # set origin for axes
# add labels
epicurve_labels +
# add theme
epicurve_theme
nice_plot
# Modify nice_plot to show only 6 breaks in the x-axis
nice_plot + scale_x_incidence(epicurve_weekly, n_breaks = 6)Now differentiating the cases by gender, using the groups = argument in the incidence command
# Create epiweek object, with counts grouped by gender
epicurve_weekly_gender <- incidence(linelist$date_of_onset,
interval = "week",
groups = linelist$gender,
na_as_group = FALSE) # Prevents missing values from being assigned their own group
# Plot the epicurve
# Note: Remove the boxes around each case as it makes gender colours hard to see! (show_cases = FALSE)
nice_plot <- plot(epicurve_weekly_gender, show_cases = FALSE, border = "black", n_breaks = nrow(epicurve_weekly_gender)) +
# add labels (defined in previous section)
epicurve_labels +
# add theme elements
epicurve_theme
nice_plotTo filter data, This version is filtered to only show data from a specific province.
# filter the dataset and pass it to the incidence() function
Zhejiang_data <- filter(linelist, province == "Zhejiang")
epicurve_Zhejiang <- incidence(Zhejiang_data$date_of_onset,
interval = "week",
groups = Zhejiang_data$gender)
# Re-sets labels, changing title to reflect subset
epicurve_labels <- labs(x = "Week",
y = "Cases",
title = "H5N7 cases by week of onset in Zhejiang",
caption = paste0("Source: Linelist data from: ", data_date, "; ", missing_onset, " are missing date of onset and not shown."))
# plot as before
plot(epicurve_Zhejiang, show_cases = TRUE, border = "grey") +
# add labels (defined in previous section)
epicurve_labels +
# add theme elements
epicurve_themeBelow are tabs on using “ggplot2” package
# Daily case counts
###################
plot_daily <- ggplot(linelist, aes(x = date_of_onset)) +
# stacked bars, bined by day (1 days)
stat_bin(binwidth = 1, position="stack")
print(plot_daily)
# Weekly case counts
###################
plot_weekly <- ggplot(linelist, aes(x = date_of_onset)) +
# stacked bars, bined by week (7 days)
stat_bin(binwidth = 7, position="stack", fill = "brown")
print(plot_weekly)# Preparation
#############
# Create epiweek variable. Factor argument automatically includes all weeks in span. Numeric shows just the week number.
linelist$epiweek <- aweek::date2week(linelist$date_of_onset, factor = TRUE, numeric = TRUE)
# Calculate maximum number of cases in an epiweek, to get the maximum y-axis height (also helps with uniformity in multiple plots)
ymax <- max(summary(factor(linelist$epiweek), maxsum = length(linelist$epiweek)))
# Weekly case counts
###################
plot_weekly <- ggplot(linelist, aes(x = date_of_onset)) +
# stacked bars, bined by week (7 days)
stat_bin(binwidth = 7, position = "stack", fill = "grey", color = "black") +
# X-axis 21-day labels
scale_x_date(
# Sets date label breaks as every 3 weeks from Monday before the first case
breaks = function(x) seq.Date(from = min(linelist$date_of_onset, na.rm = T), to = max(linelist$date_of_onset, na.rm=T), by = "1 week"),
# axis limits determined by max/min + buffer
limits = c((min(linelist$date_of_onset, na.rm = T) - 8), (max(linelist$date_of_onset, na.rm = T) + 8)),
# displays as date number, then abbreviated month (e.g. 12 Oct)
date_labels = "%d-%b",
# sets origin at (0,0)
expand = c(0,0)) +
# Y-axis breaks every 5 cases
scale_y_continuous(breaks = seq(0, ymax, 5),
limits = c(0, ymax),
expand = c(0, 0)) +
# Theme specifications (axis, text, etc.)
theme(# title
plot.title = element_text(size=20, hjust= 0, face="bold"), # title size, font, bold
# axes
axis.text.x = element_text(angle=90, vjust=0.5, hjust=1),
axis.text = element_text(size=12),
axis.title = element_text(size=14, face="bold"),
axis.line = element_line(colour="black"),
# background
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
# caption (italics, on right side)
plot.caption = element_text(hjust = 0, face = "italic")
) +
guides(fill = guide_legend(reverse = TRUE, # Orders Non-active zones at end of legend
override.aes = list(size = 0.2),
ncol = 2)) + # Number of legend columns
labs(x = "Week of illness onset",
y = "Number of cases",
subtitle = "subtitle here",
caption = paste0(nrow(linelist),
" confirmed and probable cases, reported as of ", data_date, ". ",
missing_onset, " cases missing date of onset and not shown.")) +
ggtitle("Epidemic curve")
# print
print(plot_weekly)Colored by a category
# Setup
########
# Two known classes (select colors from colorbrewer2.org)
colors_overall = c("#d95f02", #
"#1b9e77",
"#7570b3") #
# Order sex variable by reverse # of cases, so plot stacks with smallest # of cases at top
linelist$gender <- factor(linelist$gender,
levels = levels(fct_rev(fct_infreq(linelist$gender))))
# Calculates maximum yaxis height for uniformity between the two graphs
ymax <- max(summary(factor(linelist$epiweek), maxsum = length(linelist$epiweek)))
# Number missing onset_date and cannot be graphed
missing_onset <- nrow(linelist[is.na(linelist$date_of_onset),])
# PLOT - BY ONSET DATE
######################
plot_defined_cats <- ggplot(linelist, aes(x = date_of_onset, fill = gender)) +
# stacked bars, width of 7 days
stat_bin(binwidth = 7, position = "stack") +
# Colors and labels of confirmed/probable
scale_fill_manual(values = rev(colors_overall),
labels = str_to_sentence(levels(factor(linelist$gender)))) +
# X-axis scale labels (not aggregation, just the labels)
scale_x_date(# Sets date label breaks as every week
breaks = function(x) seq.Date(from = min(linelist$date_of_onset, na.rm = T), to = max(linelist$date_of_onset, na.rm = T), by = "1 week"),
limits = c((min(linelist$date_of_onset, na.rm=T)), (max(linelist$date_of_onset, na.rm = T))), # axis limits determined by max/min + buffer
date_labels = "%d-%b", # displays as date # then abbreviated month (e.g. 12 Oct)
expand = c(0, 0)) + # sets origin at (0,0)
# Y-scale in breaks, up to the ymax previously defined
scale_y_continuous(breaks = seq(0, 500, 5), limits = c(0, ymax), expand=c(0, 0)) +
# Themes for axes, titles, background, etc.
theme(plot.title = element_text(size=20, hjust=0.5, face="bold"),
axis.text = element_text(size=12),
axis.title = element_text(size=14, face="bold"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
axis.text.x = element_text(angle=90, vjust=0.5, hjust=1)) +
# Legend specifications
theme(legend.title = element_blank(),
legend.justification = c(0, 1),
legend.position = c(0.09, 0.98),
legend.background = element_blank(),
legend.text = element_text(size = 12)) +
guides(fill = guide_legend(reverse = TRUE, override.aes = list(size = 0.2))) +
# Axis and caption labels
labs(x = "Week of illness onset",
y = "Number of Cases",
caption = paste(missing_onset,"cases were missing onset date and are not included in the onset graph")) +
# Title
ggtitle("Cases by week of illness onset")
# print
print(plot_defined_cats)This code creates an epicurve with the following specifications:
# PARAMETERS
#############
# Maximum y-value for epiweek (this will be larger than necessary because of missing onset dates)
ymax <- max(table(linelist$epiweek))
# Number missing onset_date and cannot be graphed
missing_onset <- nrow(filter(linelist, is.na(date_of_onset)))
# SETUP - ACTIVE/NON-ACTIVE ZONES
#################################
# List of "active" zones with a case in the date range
active_zones <- unique(linelist$province[which(linelist$date_of_onset > (data_date - 90))])
active_zones
## [1] "Fujian" "Jiangxi" "Beijing" "Hebei" "Guangdong"
# Table of active zones and their overall number of cases (for ordering their stacked appearance)
order_table <- linelist %>%
filter(province %in% active_zones) %>%
group_by(province) %>%
summarise(cases = n())
order_table
## # A tibble: 5 x 2
## province cases
## <chr> <int>
## 1 Beijing 3
## 2 Fujian 5
## 3 Guangdong 1
## 4 Hebei 1
## 5 Jiangxi 6
# Create TRUE/FALSE variable for "active" health zones
linelist$active_zone <- ifelse(linelist$province %in% active_zones, TRUE, FALSE)
# Create list of non-active HZ names for bottom of plot
other_zone_names <- unique(sort(linelist$province[linelist$active_zone == FALSE]))
# Make variable for graph categories, including a level for "non-active" zones
linelist$graph_zone <- factor(case_when(
# Value assignments
# Non-active zones
linelist$active_zone == FALSE ~ "Non-active zones",
# All others are assigned their names, capitalized
TRUE ~ stringr::str_to_title(linelist$province)),
# Order of variable levels
levels = c(
# "Non-active zones" is first level
"Non-active zones",
# Orders active zones by their frequency in linelist, reversed, so most-affected zones are on the BOTTOM of plot
str_to_title(rev(levels(fct_infreq(as.factor(linelist$province[linelist$active_zone == TRUE])))))))
table(linelist$graph_zone, useNA = "ifany")
##
## Non-active zones Hebei Guangdong Beijing
## 120 1 1 3
## Fujian Jiangxi
## 5 6
# COLORS
########
# Number of unique values in graph_zone variable, minus 1 (for non-active, which is added later as grey (#cccccc))
colors_needed <- length(unique(linelist$graph_zone, na.rm=T)) - 1
# List of possible colors (see colorbrewer2.com, qualitative scheme)
colors_linelist = c(#"#cccccc", # first = non-active grey color
"#1b9e77", # turquoise green
"#ff7f00", # orange
"#ffff33", # yellow
"#6a3d9a", # purple
"#b15928", # brown
"#1f78b4", # blue
"#e31a1c", # red,
"#fb9a99", # pink
"#b2df8a", # light green
"#cab2d6", # light purple
"#a6cee3", # light blue
"#fdbf6f", # beige
"#33a02c" # green
)
# Reduce number of colors to only the number needed
colors_linelist <- c("#cccccc", rev(colors_linelist[1:colors_needed]))
# MAKE GRAPH
#############
plot_overall <- ggplot(linelist, aes(x = date_of_onset, fill = graph_zone)) +
# stacked bars, bined by week (7 days)
stat_bin(binwidth = 7, position = "stack") +
# Fill of bars
scale_fill_manual(values = colors_linelist,
labels = str_to_sentence(levels(factor(linelist$graph_zone)))) +
# X-axis 21-day labels
scale_x_date( # Sets date label breaks as every 3 weeks from Monday before the first case
breaks = function(x) seq.Date(from = min(linelist$date_of_onset, na.rm = T), to = max(linelist$date_of_onset, na.rm = T), by = "1 week"),
limits = c((min(linelist$date_of_onset, na.rm = T) - 8), (max(linelist$date_of_onset, na.rm = T) + 8)), # axis limits determined by max/min + buffer
date_labels = "%d-%b", # displays as date number, then abbreviated month (e.g. 12 Oct)
expand = c(0,0)) + # sets origin at (0,0)
# Y-axis breaks every 5 cases
scale_y_continuous(breaks = seq(0, ymax, 5),
limits = c(0, ymax),
expand = c(0, 0)) +
# Theme specifications (axis, text, etc.)
theme(plot.title = element_text(size = 20, hjust = 0, face = "bold"), # title size, font, bold
axis.text.x = element_text(angle=90, vjust=0.5, hjust=1),
axis.text = element_text(size=12),
axis.title = element_text(size=14, face="bold"),
axis.line = element_line(colour="black"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
plot.caption = element_text(hjust = 0, face = "italic")
) +
# Legend specifications
theme(legend.title = element_blank(), # No legend title
legend.position = c(0.20, 0.85), # placement of legend
legend.background = element_blank(), # legend background
legend.text = element_text(size=12)) + # legend text size
guides(fill = guide_legend(reverse = TRUE, # Orders Non-active zones at end of legend
override.aes = list(size = 0.2),
ncol = 2)) + # Number of legend columns
labs(x = "Week of illness onset",
y = "Number of cases",
subtitle = "Health zones with cases in the last 42 days specified by color",
caption = paste0(nrow(linelist),
" confirmed and probable cases, reported as of ", data_date, ". ",
missing_onset, " cases missing date of onset and not shown.",
"\nNon-active zones include: ", str_to_title(toString(unique(linelist$province[linelist$active_zone == FALSE]))))) +
ggtitle("Epidemic curve by active health zones")
# print
plot_overall
#SETUP
#############
# Filter to health zone of interest
zone_data <- linelist
# Number missing onset_date and cannot be graphed
missing_onset <- nrow(filter(linelist, is.na(date_of_onset)))
# Assign health area groups (individual for HAs of interest, groups others together)
linelist$graph_areas <- factor(case_when(
linelist$province == "Shanghai" ~ "Shanghai",
linelist$province == "Jiangsu" ~ "Jiangsu",
linelist$province == "Zhejiang" ~ "Zhejiang",
TRUE ~ "Other (10)"
),
# Levels part of the factor function assigns order of appearance
levels = c(
"Other (10)",
"Shanghai",
"Jiangsu",
"Zhejiang"
)
)
# checks
table(linelist$graph_areas, useNA = "ifany")
##
## Other (10) Shanghai Jiangsu Zhejiang
## 29 33 28 46
# Color assignments
colors_needed <- length(unique(linelist$graph_areas, na.rm=T)) - 1 # number of colors needed
# list of colors
colors_aire = c("#a6cee3",
"#1f78b4",
"#b2df8a",
"#33a02c",
"#fb9a99",
"#e31a1c",
"#fdbf6f",
"#ff7f00",
"#cab2d6",
"#6a3d9a",
"#ffff99",
"#b15928"
)
# Reduce number of colors to only the number needed
colors_aire <- c("#cccccc", rev(colors_aire[1:colors_needed]))
# Plot of province
#####################################
plot <- ggplot(linelist, aes(x = date_of_onset, fill = graph_areas)) +
stat_bin(binwidth = 7, position="stack") +
scale_fill_manual(values = colors_aire, labels = str_to_sentence(levels(factor(linelist$graph_areas)))) +
scale_x_date(date_breaks = "1 week", date_labels = "%d-%b", limits = c((min(linelist$date_of_onset, na.rm = T) - 8), (max(linelist$date_of_onset, na.rm = T) + 8)), expand=c(0,0)) + # I used the date onset variable here so x axes will be the same
scale_y_continuous(breaks = seq(0, 500, 5), limits = c(0, 35), expand = c(0, 0)) +
theme(plot.title = element_text(size = 20, hjust = 0.5, face = "bold"),
plot.caption = element_text(hjust = 0, face = "italic"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
axis.text = element_text(size = 12),
axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
axis.title = element_text(size = 14, face = "bold"),
legend.title = element_blank(),
legend.justification = c(0,1),
legend.position = c(0.05, 1),
legend.background = element_blank(),
legend.text = element_text(size = 12)) +
guides(fill = guide_legend(reverse = TRUE, override.aes = list(size = 0.2), ncol = 4)) +
labs(x="Week of illness onset",
y="Number of cases",
subtitle = "",
caption = paste0(nrow(zone_data), " confirmed and probable cases, as of ", data_date, ". \n", missing_onset, " cases excluded due to missing date of onset.")) +
ggtitle("Cases of influenza, by province")
plot# Define the waves
##################
# zone_data <- filter(linelist, zone_de_sante == "mabalako")
#
# zone_data$wave <- case_when(
# zone_data$date_onset >= as.Date("2018-03-01") &
# zone_data$date_onset < as.Date("2018-10-25") ~ "Wave 1",
#
# zone_data$date_onset >= as.Date("2018-10-25") &
# zone_data$date_onset < as.Date("2019-02-01") ~ "Wave 2",
#
# zone_data$date_onset >= as.Date("2019-02-01") &
# zone_data$date_onset < as.Date("2019-09-15") ~ "Wave 3",
#
# zone_data$date_onset >= as.Date("2019-09-15") ~ "Wave 4",
#
# TRUE ~ NA_character_
# )
#
# table(is.na(zone_data$date_onset))
# table(zone_data$wave, useNA = "always")
#
#
# # Color assignments
# colors_needed <- length(unique(zone_data$wave, na.rm=T)) # number of colors needed
#
# # list of colors
# colors_aire = c("#a6cee3",
# "#1f78b4",
# "#b2df8a",
# "#33a02c",
# "#fb9a99",
# "#e31a1c",
# "#fdbf6f",
# "#ff7f00",
# "#cab2d6",
# "#6a3d9a",
# "#ffff99",
# "#b15928"
# )
#
# # Reduce number of colors to only the number needed
# colors_aire <- c(rev(colors_aire[1:colors_needed]))
#
#
# # Plot of health zone colored by wave
# #####################################
# plot_Mabalako <- ggplot(zone_data, aes(x = date_onset, fill = wave)) +
#
# stat_bin(binwidth = 7, position = "stack") +
#
# scale_fill_manual(values = rev(colors_aire), labels = str_to_sentence(levels(factor(zone_data$wave)))) +
#
# scale_x_date(date_breaks = "21 days", date_labels = "%d-%b",
# limits = c((min(zone_data$date_onset, na.rm = T) - 8), (max(zone_data$date_report, na.rm = T) + 8)), expand = c(0,0)) +
#
# scale_y_continuous(breaks = seq(0, 500, 5), limits = c(0, 35), expand = c(0, 0)) +
#
# theme(text = element_text(family = "Segoe Condensed"),
# axis.text.x = element_text(angle = 90, vjust = 0.5, hjust = 1),
# axis.text = element_text(size = 12),
# axis.title = element_text(size = 14, face = "bold"),
# axis.line = element_line(colour = "black"),
#
# plot.title = element_text(size = 20, hjust = 0.5, face = "bold"),
# plot.caption = element_text(hjust = 0, face = "italic"),
#
# panel.grid.major = element_blank(),
# panel.grid.minor = element_blank(),
# panel.background = element_blank(),
#
# legend.title = element_blank(),
# legend.justification = c(0,1),
# legend.position = c(0.75, 0.98),
# legend.background = element_blank(),
# legend.text = element_text(size=12)) +
#
# guides(fill = guide_legend(reverse = TRUE, override.aes = list(size = 0.2), ncol = 1)) +
#
# labs(x="Week of illness onset",
# y="Number of cases",
# subtitle = "",
# caption = paste0(nrow(zone_data), " confirmed and probable cases, as of ", data_date, ". \n", missing_onset, " cases excluded due to missing date of onset and 16 excluded due to uncertain health zone of report.")) +
#
# ggtitle("Four waves of EVD in Mabalako health zone")
#
# plot_Mabalako
#
#
# # Produce table describing each wave
# ####################################
# table <- zone_data %>%
# select("aire_de_sante", "wave", "community_death", "date_onset", "cte_date", "epicasedef", "community_death", "contact_registered", "contact_surveilled") %>%
# group_by(wave) %>%
# summarise(first_onset = min(date_onset, na.rm = T),
# last_admission = max(cte_date, na.rm = T),
# n = n(),
# confirmed = sum(epicasedef == "confirmed"),
# community_deaths = paste0(sum(community_death == 1),
# " (", round(100*sum(community_death == 1)/confirmed),"%)"),
# reg_contacts = paste0(sum(contact_registered == "yes"),
# " (", round(100*sum(contact_registered == "yes")/confirmed),"%)"),
# surv_contacts = paste0(sum(contact_surveilled == "yes"),
# " (", round(100*sum(contact_surveilled == "yes")/confirmed),"%)"),
# top = paste(toupper(names(sort(table(aire_de_sante),decreasing=TRUE)[1:3])), collapse=", ",
# round(100*(sort(table(aire_de_sante),decreasing=TRUE)[1:3]/confirmed)), "%"),
# health_areas = paste(toupper(unique(aire_de_sante)), collapse=', ')
# )
#
# kable(table)Often you do not have linelist data, but instead daily case counts from facilities, districts, etc. You can plot these in an epidemiological curve, but the code will be slightly different.
This section will utilize the counts_data dataset that was imported earlier, in the data preparation section.
Note: The incidence package does not support aggregate data
As before, we must ensure date variables are correctly classified.
# Create epiweek variable
# aweek weeks are also stored as dates, facilitating better display manipulation
count_data$epiweek <- aweek::date2week(count_data$Date, # use the Date variable
week_start = "Monday", # epiweek begins on Monday
floor_day = TRUE, # only display year and week #
factor = TRUE) # expand to include all possible weeksggplot(data = count_data, aes(x = as.Date(epiweek), y = Cases, group = District, fill = District))+
geom_bar(stat = "identity")+
# LABELS for x-axis
scale_x_date(date_breaks = "1 month", # displays by month
date_labels = '%b%d\n%Y')+ #labeled by month with year below
# Choose color palette (uses RColorBrewer package)
scale_fill_brewer(palette = "Pastel1")+
# Theme specifications (axis, text, etc.)
theme(
# title
plot.title = element_text(size=20, hjust= 0, face="bold"), # title size, font, bold
# axes
axis.text.x = element_text(angle=0, vjust=0.5, hjust=1),
axis.text = element_text(size=12),
axis.title = element_text(size=14, face="bold"),
axis.line = element_line(colour="black"),
# background
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
# caption (italics, on right side)
plot.caption = element_text(hjust = 0, face = "italic"))+
# labels
labs(x = "Week of report",
y = "Number of cases",
subtitle = "Cases aggregated by week and shown by district",
caption = "Data source: XXXXX")+
ggtitle("Epidemic curve of disease X in fictional location")Although there are fierce discussions about the validity of this within the data visualization community, many supervisors want to see an epicurve or similar chart with a percent overlaid with a second axis.
In ggplot it is very difficult to do this, except for the case where you are showing a line reflecting the proportion of a category shown in the bars below.
This uses the linelist dataset
TODO not complete yet
library(reshape2)
# group the data by week, summarize counts by group (gender)
linelist_week <- linelist %>%
mutate(onset_epiweek = aweek::date2week(date_of_onset, floor_day = TRUE, factor = TRUE)) %>%
group_by(onset_epiweek) %>%
summarize(num_male = sum(gender == "m"),
num_female = sum(gender == "f"),
pct_male = round(100*(num_male / n())),
med_age = median(as.numeric(age), na.rm=T)
)
# remove pct and melt
linelist_week_melted <- linelist_week %>%
select(-c("pct_male", "med_age")) %>%
melt(id.vars = c("onset_epiweek"))
# merge together (multiple of the same values in week will attach to melted)
linelist_week_melted <- merge(linelist_week_melted,
linelist_week,
by = "onset_epiweek")
second_axis <- ggplot(linelist_week_melted,
aes(x = as.Date(onset_epiweek),
y = value, group = variable,
fill = variable)) +
# bars
geom_bar(stat = "identity")+
# Colors and labels of confirmed/probable
scale_fill_manual(values = c("blue", "red"),
labels = str_to_sentence(levels(factor(linelist_week_melted$variable)))) +
geom_line(mapping = aes(y = pct_male, color = "% male"), size = 0.5) +
scale_color_manual(values = "black")+
scale_y_continuous(sec.axis = sec_axis(~(./sum(linelist_week_melted$value, na.rm = T)*100), name = "name here", breaks = seq(0, 100, 20)))+
# X-axis scale labels (not aggregation, just the labels)
scale_x_date(# Sets date label breaks as every week
breaks = function(x) seq.Date(from = min(linelist$date_of_onset, na.rm = T), to = max(linelist$date_of_onset, na.rm = T), by = "1 week"),
limits = c((min(linelist$date_of_onset, na.rm=T)), (max(linelist$date_of_onset, na.rm = T))), # axis limits determined by max/min + buffer
date_labels = "%d-%b", # displays as date # then abbreviated month (e.g. 12 Oct)
expand = c(0, 0)) + # sets origin at (0,0)
# Y-scale in breaks, up to the ymax previously defined
scale_y_continuous(breaks = seq(0, 500, 5), limits = c(0, ymax), expand=c(0, 0)) +
# Themes for axes, titles, background, etc.
theme(plot.title = element_text(size=20, hjust=0.5, face="bold"),
axis.text = element_text(size=12),
axis.title = element_text(size=14, face="bold"),
panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
axis.line = element_line(colour = "black"),
axis.text.x = element_text(angle=90, vjust=0.5, hjust=1)) +
# Legend specifications
theme(legend.title = element_blank(),
legend.justification = c(0, 1),
legend.position = c(0.09, 0.98),
legend.background = element_blank(),
legend.text = element_text(size = 12)) +
guides(fill = guide_legend(reverse = TRUE, override.aes = list(size = 0.2))) +
# Axis and caption labels
labs(x = "Week of illness onset",
y = "Number of Cases",
caption = paste(missing_onset,"cases were missing onset date and are not included in the onset graph")) +
# Title
ggtitle("Cases by week of illness onset")
second_axis
# print
print(plot_defined_cats)Here are links to good online resources
Other resources
This analysis plots the frequency of different combinations of values/responses. In this example, we plot the frequency of different symptom combinations.
This analysis is often called:
Multiple response analysis Sets analysis Combinations analysis
The first method shown uses the package ggupset. More information can be found online or offline in the package documentation in your RStudio Help tab.
This linelist includes five “yes/no” variables on reported symptoms. We will need to transform these variables a bit to use the ggupset package to make our plot.
View the data (scroll to the right to see the symptoms variables)
We convert the “yes” and “no the the actual symptom name. If”no", we set the value as blank.
# create column with the symptoms named, separated by semicolons
linelist_sym_1 <- linelist_sym %>%
# convert the "yes" and "no" values into the symptom name itself
mutate(fever = case_when(fever == "yes" ~ "fever", # if old value is "yes", new value is "fever"
TRUE ~ NA_character_), # if old value is anything other than "yes", the new value is NA
chills = case_when(chills == "yes" ~ "chills",
TRUE ~ NA_character_),
cough = case_when(cough == "yes" ~ "cough",
TRUE ~ NA_character_),
aches = case_when(aches == "yes" ~ "aches",
TRUE ~ NA_character_),
shortness_of_breath = case_when(shortness_of_breath == "yes" ~ "shortness_of_breath",
TRUE ~ NA_character_))Now we make two final variables:
1. Pasting together all the symptoms of the patient (character variable)
2. Convert the above to class list, so it can be accepted by ggupset to make the plot
linelist_sym_1 <- linelist_sym_1 %>%
mutate(
# combine the variables into one, using paste() with a semicolon separating any values
all_symptoms = paste(fever, chills, cough, aches, shortness_of_breath, sep = "; "),
# make a copy of all_symptoms variable, but of class "list" (which is required to use ggupset() in next step)
all_symptoms_list = as.list(strsplit(all_symptoms, "; "))
)View the new data. Note the two columns at the end - the pasted combined values, and the list
Load required package to make the plot (ggupset)
Create the plot:
ggplot(linelist_sym_1,
aes(x=all_symptoms_list)) +
geom_bar() +
scale_x_upset(reverse = FALSE,
n_intersections = 10,
sets = c("fever", "chills", "cough", "aches", "shortness_of_breath")
)+
labs(title = "Signs & symptoms",
subtitle = "10 most frequent combinations of signs and symptoms",
caption = "Caption here.",
x = "Symptom combination",
y = "Frequency in dataset")The UpSetR package allows more customization, but it more difficult to execute:
https://github.com/hms-dbmi/UpSetR read this https://gehlenborglab.shinyapps.io/upsetr/ Shiny App version - you can upload your own data https://cran.r-project.org/web/packages/UpSetR/UpSetR.pdf documentation - difficult to interpret
Convert symptoms variables to 1/0.
# Make using upSetR
linelist_sym_2 <- linelist_sym %>%
# convert the "yes" and "no" values into the symptom name itself
mutate(fever = case_when(fever == "yes" ~ 1, # if old value is "yes", new value is "fever"
TRUE ~ 0), # if old value is anything other than "yes", the new value is NA
chills = case_when(chills == "yes" ~ 1,
TRUE ~ 0),
cough = case_when(cough == "yes" ~ 1,
TRUE ~ 0),
aches = case_when(aches == "yes" ~ 1,
TRUE ~ 0),
shortness_of_breath = case_when(shortness_of_breath == "yes" ~ 1,
TRUE ~ 0))Now make the plot, using only the symptom variables. Must designate which “sets” to compare (the names of the symptom variables).
Alternatively use nsets = and order.by = "freq" to only show the top X combinations.
# Make the plot
UpSetR::upset(
select(linelist_sym_2, fever, chills, cough, aches, shortness_of_breath),
sets = c("fever", "chills", "cough", "aches", "shortness_of_breath"),
order.by = "freq",
sets.bar.color = c("blue", "red", "yellow", "darkgreen", "orange"), # optional colors
empty.intersections = "on",
# nsets = 3,
number.angles = 0,
point.size = 3.5,
line.size = 2,
mainbar.y.label = "Symptoms Combinations",
sets.x.label = "Patients with Symptom")Also, 3d?
via ggplot and via R4Epis methods
ggplot()List of different diagrams in these tabs…
E.g. EVD patient “pathways” to outcome (via clinic or not, etc.)
HIV care continuum datasets? PreP datasets?
Sankey plots - show transitions among cohort over time, interrelatedness of groups Liza Coyer TODO
Or papers in meta-analysis
E.g. border closures during COVID
Why How When etc.
Many audiences and reasons for the tables…
knitr::kable DT
For publication
quickly changing the denominator (per 100,000, etc.)
Embed Rmarkdown cheatsheet Table issues HTML candies Making tables, cheatsheet contained in HTML handbook?
Endemic corridor analysis Detecting spikes in syndromic/routine surveillance
{.tabset .tabset-fade .tabset-pills} stringr, gsub,
antijoins, etc dplyr
rowmatcher Other methods
Identifying and getting rid of duplicates
Combining the values of multiple records into one record
.Rprofile
overview
here package
How to quickly exclude or include in tables, charts, analysis, exclude whole line if any is missing, all missing, etc.
Embed ggplot cheatsheet
gghighlight
ggrepel
Cowplot Complicated method (% 100 * …)
Tidymodels
Liza Coyer TODO this, logitudinal
Flexdashboard
Troubleshooting tips, common errors, etc.
plotly
If there is sufficient file space, this will be a gallery of visualizations, tables, etc. with internal links to the relevant Handbook sections.